En este apunte mostramos cómo utilizar JAGS para inferir los valores de posición $\mu$ y de dispersión $\sigma$ de una distribución Normal.
Comenzaremos simulando ciertas observaciones Normales con R a partir de valores conocidos de $\mu$ y $\sigma$. Después utilizaremos esas observaciones en JAGS para evaluar si dicho programa puede recuperar los valores paramétricos que, sabemos, "generaron" las observaciones. Estos dos pasos se conocen como la técnica de "recuperación paramétrica", y se utiliza para evaluar el funcionamiento de algún modelo bajo algún paradigma de inferencia (no necesariamente Bayesiano).
En la célula siguiente le pedimos a R que genere observaciones provenientes de una distribución Normal con parámetros $\mu=53$ y $\sigma=10$. Noten que en esta simulación el parámetro de dispersión que R espera (en las líneas 4 y 6 de la célula) está en unidades de desviación estándar. Después de simular las observaciones y graficar, incluimos algunos descriptivos de nuevo arreglo z; junto con el histograma, estos primeros pasos son muy importantes para comenzar a estudiar las características de cualquier conjunto de datos
rm(list=ls())
true_mu <-53
true_sigma <- 10 # Desviación Estándar
z_plot <- seq(0,100,length.out=100)
pdf_z <- dnorm(z_plot,mean=true_mu,sd=true_sigma) # Al calcular la densidad normal en R, la función espera DESVIACIÓN ESTÁNDAR
set.seed(923) # Esta instrucción permite reproducir el proceso pseudo-aleatorio que genera observaciones al "azar" desde cierto modelo
simul_z <- rnorm(n=50,mean=true_mu,sd=true_sigma) # Al generar observaciones normales en R, la función espera DESVIACIÓN ESTÁNDAR
# Histograma de las obs. simuladas
col_data <- '#fad139'
brks <- seq(0,100,3)
options(repr.plot.width = 10, repr.plot.height = 7)
hist(simul_z,xlim=c(0,100),breaks=brks,freq=F,border=F,col=paste(col_data,'88',sep=''),ann=F,axes=F)
points(simul_z,rep(0,length(simul_z)),pch=21,bg='#000000',col=col_data,cex=3,lwd=3)
lines(z_plot,pdf_z,lwd=5,col=col_data)
axis(1,cex.axis=2)
mtext('z',1,cex=3,line=3)
mtext('Observaciones',3,cex=2)
summary(simul_z)
Min. 1st Qu. Median Mean 3rd Qu. Max. 34.18 47.53 51.36 52.38 57.29 72.33
La idea de simular observaciones es que ahora conocemos exactamente qué modelo generó los valores en simul_z, y no sólo eso, sino que sabemos cuáles eran los valores "verdaderos" de $\mu$ y de $\sigma$ "responsables" de esas observaciones.
Ahora intentaremos recuperar dichos valores paramétricos programando el mismo modelo (Normal) en JAGS y analizando las distribuciones posteriores correspondientes. En concreto, utilizaremos simul_z como si se trataran de mediciones recolectadas en el mundo natural e intentaremos entenderlas utlizando el modelo:

En otras palabras, comenzaremos suponiendo que el centro de la distribución $\mu$ (también llamado "media") puede estar en cualquier posición entre $-100$ y $100$, y que el parámetro $\sigma$ que representa desviación estándar puede ser cualquier valor entre $0$ y $50$.
Para programar este modelo en JAGS seguiremos la misma estructura que en modelos anteriores, aunque poniendo especial atención a la hora de expresar el parámetro de dispersión en el modelo: a diferencia de R, JAGS no parametriza a la distribución Normal en términos de desviación estándar, sino en términos de precisión $\tau$, en donde la relación entre ambas es:
o bien:
$$\begin{align} \sigma^2\tau&=\frac{\sigma^2}{\sigma^2}\\\\ \sigma^2\tau&=1\\\\ \sigma^2&=\frac{1}{\tau}\\\\ \sqrt{\sigma^2}&=\sqrt{\frac{1}{\tau}}\\\\ \sigma&=\frac{\sqrt{1}}{\sqrt{\tau}}\\\\ \sigma&=\frac{1}{\sqrt{\tau}} \end{align}$$es decir, en un modelo normal lo que se conoce como precisión es igual al recíproco de la varianza, o bien al recíproco de la desviación estándar al cuadrado. Las tres etiquetas, o los tres parámetros ($\tau$, $\sigma$, y $\sigma^2$) sirven para medir la dispersión de la curva, aunque cada uno opera en escalas diferentes. Lo importante es que la relación entre ellos es precisa y podemos programarla para extraer e interpretar el parámetro de dispersión en la escala deseada:
# Unobserved
unobserved <- c('mu_post',
'mu_prior',
'tau_post',
'tau_prior',
'sd_post',
'sd_prior',
'z_post', # Post. Predictiva (a.k.a. "Postdicción")
'z_prior') # Prior Predictiva
# Observed
z <- simul_z
n_obs <- length(z) # Cuando se trabaja con más de una observación es necesario
# hacer explícito (para JAGS) la cantidad total de observaciones .
observed <- list('z','n_obs') # Noten que incluimos el número de observaciones
# en esta lista.
# Model
write('
model{
# Unobserved (priors)
mu_post~dunif(-100,100)
mu_prior~dunif(-100,100)
sd_post~dunif(0,60) # Expresamos incertidumbre inicial en escala de "desv. est."...
sd_prior~dunif(0,60)
tau_post <- 1/sd_post^2 # ...y después traducimos a escala de "precisión"
tau_prior <- 1/sd_prior^2
# Observed (likelihood)
for(j in 1:n_obs){ # La notación de platos se traduce en un ciclo FOR en el modelo
z[j]~dnorm(mu_post,tau_post) # JAGS espera PRECISIÓN en la Normal (i.e., DEVOLVERÁ precisión
# sin importar el *nombre* del segundo argumento en esta función)
}
# Predictive
z_prior~dnorm(mu_prior,tau_prior)
z_post~dnorm(mu_post,tau_post)
}
','gaussian_model_jags.bug')
# Inference
library('R2jags')
set.seed(95)
bayes <- jags(data = observed,
parameters.to.save = unobserved,
model.file = 'gaussian_model_jags.bug',
n.iter=3000,
n.chains=4,
n.thin=3,
n.burnin=0)
nodos <- bayes$BUGSoutput$sims.list
Compiling model graph Resolving undeclared variables Allocating nodes Graph information: Observed stochastic nodes: 50 Unobserved stochastic nodes: 6 Total graph size: 67 Initializing model
En otras palabras, para escribir el modelo original en JAGS es necesario (y será necesario siempre que trabajemos con modelos Normales) programar esta versión ligeramente modificada:
$$\begin{align} \mu&\sim Uniform(-100,100)\\ \sigma&\sim Uniform(0,60)\\ \tau&\leftarrow\frac{1}{\sigma^2}\\ z_j&\overset{iid}{\sim}Gaussian(\mu,\tau) \end{align}$$
bayes$BUGSoutput$summary
| mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | Rhat | n.eff | |
|---|---|---|---|---|---|---|---|---|---|
| deviance | 356.99176681 | 2.07136328 | 3.549771e+02 | 3.555253e+02 | 356.352903663 | 3.578048e+02 | 362.65085579 | 1.001321 | 3000 |
| mu_post | 52.40581578 | 1.23388370 | 4.993261e+01 | 5.159479e+01 | 52.421931477 | 5.323664e+01 | 54.79885243 | 1.002226 | 1500 |
| mu_prior | -0.90712295 | 58.48263283 | -9.568755e+01 | -5.204934e+01 | -0.419440182 | 4.914262e+01 | 95.27402921 | 1.001801 | 1900 |
| sd_post | 8.71203666 | 0.90261550 | 7.165102e+00 | 8.077755e+00 | 8.625855272 | 9.265346e+00 | 10.70420007 | 1.000776 | 4000 |
| sd_prior | 30.14921657 | 17.47146232 | 1.649450e+00 | 1.503261e+01 | 30.141797618 | 4.542286e+01 | 58.42027384 | 1.000769 | 4000 |
| tau_post | 0.01359167 | 0.00274725 | 8.727534e-03 | 1.164868e-02 | 0.013439888 | 1.532564e-02 | 0.01947849 | 1.000776 | 4000 |
| tau_prior | 0.78588661 | 19.46098092 | 2.930035e-04 | 4.846755e-04 | 0.001100682 | 4.425182e-03 | 0.36755432 | 1.000769 | 4000 |
| z_post | 52.46051563 | 8.92099768 | 3.497255e+01 | 4.640622e+01 | 52.513684795 | 5.819371e+01 | 70.06369460 | 1.001504 | 2700 |
| z_prior | -0.49500746 | 67.38175900 | -1.188988e+02 | -5.408884e+01 | -0.939428367 | 5.093649e+01 | 126.70906921 | 1.001294 | 3100 |
options(repr.plot.width = 14, repr.plot.height = 3)
traceplot(bayes)
Los estadísticos $\hat{R}$ y n.eff, y la apariencia visual de las cadenas, sugieren convergencia confiable. Podemos pasar a examinar los resultados de JAGS para decidir si el algoritmo logró "adivinar" los valores arbitrarios de $\mu$ y $\sigma$ con los que generamos las observaciones $z_j$.
El primer paso para examinar el resultado de JAGS siempre consiste en comparar la distribución postdictiva del modelo contra las mediciones de las variables observables. Si estas distribuciones se parecen, podemos continuar interpretando el resto de parámetros, pero si no se parecen y difieren radicalmente no es conveniente continuar en tanto que los parámetros restantes no tendrán sentido respecto de las observaciones porque el modelo que los relaciona es incorrecto.
col_prior <- '#41bbc5'
col_post <- '#1e438d'
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Zoom izquierda
brks <- seq(-300,300,3)
hist(nodos$z_post,breaks=brks,xlim=c(-150,150),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$z_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'88',sep=''))
axis(1,cex.axis=1.7)
mtext('z (predicciones)',1,cex=2,line=3)
legend(-150,max(h_pst$density),legend=c('prior predictiva','postdicción'),pch=15,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.3,pt.cex=2,box.lty='blank',bg='#00000000')
# Zoom izquierda
brks <- seq(0,100,3)
hist(simul_z,xlim=c(0,100),freq=F,breaks=brks,border=F,col=paste(col_data,'cc',sep=''),ann=F,axes=F)->h_pst
points(simul_z,rep(0,length(simul_z)),pch=21,bg='#000000',col=col_data,cex=2,lwd=3)
brks <- seq(-300,300,3)
hist(nodos$z_post,breaks=brks,xlim=c(-120,120),freq = F,add=T,border=F,col=paste(col_post,'66',sep=''))
axis(1,cex.axis=1.7)
mtext('z (observaciones)',1,cex=2,line=3)
legend(0,max(h_pst$density),legend=c('observaciones','postdicción'),pch=15,
col=c(paste(col_data,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.3,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre z',3,cex=2,outer=T,line=-2)
La distribución postdictiva del modelo se parece a la distribución de las "observaciones" $z_j$. Por lo tanto, podemos continuar examinando el resto de parámetros del modelo.
Comenzaremos examinando $\mu$, el parámetro de posición de la curva normal. Presentamos la misma gráfica con dos "zooms" diferentes para resaltar los detalles de las conclusiones posteriores. Aparte, en la gráfica de la derecha hemos incluido el "valor verdadero" de $\mu$ con el que generamos el arreglo $z_j$ para argumentar que Bayes/JAGS recupera razonablemente bien el valor de $\mu$ objetivo.
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Zoom izquierda
brks <- seq(-100,100,2)
hist(nodos$mu_post,breaks=brks,xlim=c(-120,120),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$mu_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
axis(1,cex.axis=1.7,at=seq(-120,120,30))
mtext('\u03bc',1,cex=3,line=3)
legend(-150,max(h_pst$density),legend=c('prior','posterior'),pch=15,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.4,pt.cex=2,box.lty='blank',bg='#00000000')
# Zoom derecha
brks <- seq(-100,100,.25)
hist(nodos$mu_post,breaks=brks,xlim=c(40,60),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$mu_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
abline(v=true_mu,lty='dotted',lwd=5,col='#ee0000')
axis(1,cex.axis=1.7)
mtext('\u03bc',1,cex=3,line=3)
legend(40,max(h_pst$density),legend=c('prior','posterior','\u03bc "verdadera"'),pch=c(15,15,NA),
lty=c(NA,NA,'dotted'),lwd=4,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep=''),
'#ee0000'),cex=1.5,seg.len=.5,
x.intersp=.2,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre \u03bc',3,cex=2,outer=T,line=-2)
Continuando con el parámetro de dispersión en escala de desviación estándar ($\sigma$), las gráficas siguientes permiten comparar la incertidumbre prior y posterior sobre este parámetro. Nuevamente, el "valor verdadero" de $\sigma$ se encuentra en una zona de alta densidad posterior en el resultado de Bayes/JAGS.
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Zoom izquierda
brks <- seq(-100,100,1)
hist(nodos$sd_post,breaks=brks,xlim=c(-20,80),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$sd_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
axis(1,cex.axis=1.7)
mtext('\u03c3',1,cex=3,line=3)
legend(20,max(h_pst$density),legend=c('prior','posterior'),pch=15,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.4,pt.cex=2,box.lty='blank',bg='#00000000')
# Zoom derecha
brks <- seq(-100,100,.25)
hist(nodos$sd_post,breaks=brks,xlim=c(0,20),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$sd_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
abline(v=true_sigma,lty='dotted',lwd=5,col='#ee0000')
axis(1,cex.axis=1.7)
mtext('\u03c3',1,cex=3,line=3)
legend(13,max(h_pst$density),legend=c('prior','posterior','\u03c3 "verdadera"'),pch=c(15,15,NA),
lty=c(NA,NA,'dotted'),lwd=4,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep=''),
'#ee0000'),cex=1.5,seg.len=.5,
x.intersp=.2,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre \u03c3',3,cex=2,outer=T,line=-2)
También exploraremos el parámetro de dispersión en escala de precisión ($\tau$) para enfatizar las diferencias entre esta escala y la escala de desviación estándar. Noten los valores del eje soporte y compárenlos con los valores de $\sigma$. Aparte, es importante resaltar que el suponer incertidumbre a priori uniforme sobre la desviación estándar no implica la misma forma de la incertidumbre sobre $\tau$. En este sentido, es muy importante tener presente sobre cuál de los dos parámetros expresamos nuestra incertidumbre inicial, y recordar que JAGS siempre devolverá resultados en escala de precisión. Como usuarios, es nuestro trabajo traducir dichos resultados a escala de desviación estándar.
Este ejemplo es otra forma de comprobar que Bayes/JAGS rastrea el "valor verdadero" de precisión: noten que dicho valor puede computarse a partir del valor verdadero que elegimos para la desviación estándar, utilizando la transformación en la línea 7 de la célula siguiente:
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Zoom derecha
brks <- seq(0,1000,.00075)
hist(nodos$tau_prior,breaks=brks,xlim=c(0,0.05),freq = F,border=F,col=paste(col_prior,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$tau_post,breaks=brks,add=T,freq = F,border=F,col=paste(col_post,'88',sep=''))
abline(v=1/true_sigma^2,lty='dotted',lwd=5,col='#ee0000')
axis(1,cex.axis=1.7)
mtext('\u03c4',1,cex=3,line=3)
legend(.025,max(h_pst$density),
legend=c('prior',
'posterior',
'\u03c4 "verdadera"'),
pch=c(15,
15,
NA),
lty=c(NA,
NA,
'dotted'),
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep=''),
'#ee0000'),
lwd=4,cex=1.5,seg.len=.5,
x.intersp=.2,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre \u03c4',3,cex=2,line=1)
La base de observaciones dataset_01.csv contiene registros de una variable.
R, haz un histograma de las observaciones, y extrae los descriptivos básicos.Suponiendo que las observaciones provienen de un modelo Normal, infiere los parámetros $\mu$ y $\sigma$ correspondientes.
a. Ajusta el modelo de JAGS si es necesario, sobre todo respecto de los límites apropiados para los no-observables.
b. Evalúa el resultado del algoritmo.
c. Evalúa el desempeño del modelo, comparando la distribución postdictiva contra las observaciones originales.
d. Interpreta las distribuciones posteriores sobre $\mu$ y $\sigma$.
rm(list=ls())
datos <- read.csv('~/Teaching/IAD-3/Notebooks/dataset_01.csv')
options(repr.plot.width = 10, repr.plot.height = 7)
col_data <- '#fd5917'
brks <- seq(0,300,3)
hist(datos$x,xlim=c(0,250),freq=F,breaks=brks,border=F,col=paste(col_data,'cc',sep=''),ann=F,axes=F)->h_pst
axis(1,cex.axis=1.7)
mtext('x (observaciones)',1,cex=2,line=3)
summary(datos$x)
Min. 1st Qu. Median Mean 3rd Qu. Max. 100.2 134.7 147.2 148.4 160.7 215.1
a. Los cambios respecto del modelo anterior se encuentran en las líneas 8, 9, 12, 13, 15, 23, 24, 32, 37, y 38 de la célula siguiente.
# Unobserved
unobserved <- c('mu_post',
'mu_prior',
'tau_post',
'tau_prior',
'sd_post',
'sd_prior',
'x_post', # Post. Predictiva (a.k.a. "Postdicción")
'x_prior') # Prior Predictiva
# Observed
x <- datos$x
n_obs <- length(x) # Cuando se trabaja con más de una observación es necesario
# hacer explícito (para JAGS) la cantidad total de observaciones .
observed <- list('x','n_obs') # Noten que incluimos el número de observaciones
# en esta lista.
# Model
write('
model{
# Unobserved (priors)
mu_post~dunif(0,200)
mu_prior~dunif(0,200)
sd_post~dunif(0,60) # Expresamos incertidumbre inicial en escala de "desv. est."...
sd_prior~dunif(0,60)
tau_post <- 1/sd_post^2 # ...y después traducimos a escala de "precisión"
tau_prior <- 1/sd_prior^2
# Observed (likelihood)
for(j in 1:n_obs){ # La notación de platos se traduce en un ciclo FOR en el modelo
x[j]~dnorm(mu_post,tau_post) # JAGS espera PRECISIÓN en la Normal (i.e., DEVOLVERÁ precisión
# sin importar el *nombre* del segundo argumento en esta función)
}
# Predictive
x_prior~dnorm(mu_prior,tau_prior)
x_post~dnorm(mu_post,tau_post)
}
','gaussian_model_jags.bug')
# Inference
library('R2jags')
set.seed(94)
bayes <- jags(data = observed,
parameters.to.save = unobserved,
model.file = 'gaussian_model_jags.bug',
n.iter=3000,
n.chains=4,
n.thin=3,
n.burnin=0)
nodos <- bayes$BUGSoutput$sims.list
summary(nodos$tau_prior)
Compiling model graph Resolving undeclared variables Allocating nodes Graph information: Observed stochastic nodes: 250 Unobserved stochastic nodes: 6 Total graph size: 266 Initializing model
V1 Min. : 0.00028 1st Qu.: 0.00050 Median : 0.00114 Mean : 0.21677 3rd Qu.: 0.00479 Max. :169.45781
b. Evaluando el resultado del algoritmo:
bayes$BUGSoutput$summary
# Nota: el valor de n.eff debe ser lo más cercano posible a: ((n.iter-n.burnin)/n.thin)*n.chains
| mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | Rhat | n.eff | |
|---|---|---|---|---|---|---|---|---|---|
| deviance | 2.204529e+03 | 1.965840e+00 | 2.202593e+03 | 2.203112e+03 | 2.203944e+03 | 2.205287e+03 | 2.209965e+03 | 1.001196 | 4000 |
| mu_post | 1.483786e+02 | 1.244525e+00 | 1.459169e+02 | 1.475470e+02 | 1.483994e+02 | 1.492150e+02 | 1.508096e+02 | 1.000853 | 4000 |
| mu_prior | 1.000059e+02 | 5.763254e+01 | 5.396310e+00 | 4.998996e+01 | 9.979163e+01 | 1.488707e+02 | 1.944555e+02 | 1.001174 | 3700 |
| sd_post | 1.996284e+01 | 8.998537e-01 | 1.831325e+01 | 1.933867e+01 | 1.993952e+01 | 2.054979e+01 | 2.186648e+01 | 1.001195 | 3600 |
| sd_prior | 2.973436e+01 | 1.733957e+01 | 1.509392e+00 | 1.444615e+01 | 2.960035e+01 | 4.476397e+01 | 5.835258e+01 | 1.000569 | 4000 |
| tau_post | 2.524565e-03 | 2.266949e-04 | 2.091425e-03 | 2.368019e-03 | 2.515188e-03 | 2.673910e-03 | 2.981735e-03 | 1.001195 | 3600 |
| tau_prior | 2.167730e-01 | 3.463047e+00 | 2.936837e-04 | 4.990487e-04 | 1.141317e-03 | 4.791766e-03 | 4.389308e-01 | 1.000569 | 4000 |
| x_post | 1.481260e+02 | 1.977677e+01 | 1.092017e+02 | 1.344183e+02 | 1.484037e+02 | 1.615518e+02 | 1.865271e+02 | 1.000679 | 4000 |
| x_prior | 1.007186e+02 | 6.734347e+01 | -2.064451e+01 | 4.936124e+01 | 1.007064e+02 | 1.516726e+02 | 2.235102e+02 | 1.002120 | 1500 |
options(repr.plot.width = 10, repr.plot.height = 3)
traceplot(bayes)
c. Evaluando el desempeño del modelo:
col_data <- '#fd5917'
col_prior <- '#bd58db'
col_post <- '#621da6'
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Zoom izquierda
brks <- seq(-350,350,3)
hist(nodos$x_post,breaks=brks,xlim=c(0,250),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$x_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'88',sep=''))
axis(1,cex.axis=1.7)
mtext('x (predicciones)',1,cex=2,line=3)
legend(0,max(h_pst$density),legend=c('prior predictiva','postdicción'),pch=15,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.3,pt.cex=2,box.lty='blank',bg='#00000000')
# Zoom izquierda
brks <- seq(0,300,3)
hist(datos$x,xlim=c(0,250),freq=F,breaks=brks,border=F,col=paste(col_data,'cc',sep=''),ann=F,axes=F)->h_pst
points(datos$x,rep(0,length(datos$x)),pch=21,bg='#000000',col=col_data,cex=2,lwd=3)
brks <- seq(-300,300,3)
hist(nodos$x_post,breaks=brks,xlim=c(-120,120),freq = F,add=T,border=F,col=paste(col_post,'66',sep=''))
axis(1,cex.axis=1.7)
mtext('x (observaciones)',1,cex=2,line=3)
legend(0,max(h_pst$density),legend=c('observaciones','postdicción'),pch=15,
col=c(paste(col_data,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,
x.intersp=.3,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre x',3,cex=2,outer=T,line=-2)
c. Interpretando la incertidumbre sobre las variables no-observables en el modelo:
options(repr.plot.width = 18, repr.plot.height = 7)
layout(matrix(1:2,ncol=2))
# Posterior sobre MU
#brks <- seq(0,200,.25)
brks <- seq(0,200,1)
hist(nodos$mu_post,breaks=brks,xlim=c(100,200),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$mu_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
axis(1,cex.axis=1.7)
mtext('\u03bc',1,cex=3,line=3)
legend(100,max(h_pst$density),legend=c('prior','posterior'),pch=c(15,15,NA),
lwd=NA,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,seg.len=.5,
x.intersp=.2,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre \u03bc',3,cex=2,outer=F,line=0)
# Posterior sobre SIGMA (desviación estándar)
brks <- seq(0,100,.6)
hist(nodos$sd_post,breaks=brks,xlim=c(0,70),freq = F,border=F,col=paste(col_post,'88',sep=''),ann=F,axes=F)->h_pst
hist(nodos$sd_prior,breaks=brks,add=T,freq = F,border=F,col=paste(col_prior,'aa',sep=''))
axis(1,cex.axis=1.7)
mtext('\u03c3',1,cex=3,line=3)
legend(35,max(h_pst$density),legend=c('prior','posterior'),pch=c(15,15),
lwd=NA,
col=c(paste(col_prior,'bb',sep=''),
paste(col_post,'bb',sep='')),cex=1.5,seg.len=.5,
x.intersp=.2,pt.cex=2,box.lty='blank',bg='#00000000')
mtext('Incertidumbre sobre \u03c3',3,cex=2,outer=F,line=0)
A continuación presentamos algunas notas sobre la función hist() en R base, con énfasis sobre cómo controlar el ancho de barras y cómo extraer la información del histograma para graficarla en otros formatos.
Para controlar los ejes y las etiquetas recomendamos examinar los argumentos ann=F y axes=F en las gráficas del apunte, así como las funciones axis() y mtext() en dichas células.
# El argumento "breaks=" permite controlar el ancho de barras con precisión:
brks <- seq(from=0.001,to=0.004,by = 0.00005) # Definimos dónde están los cortes de barra
hist(nodos$tau_post,xlim=c(0,.005),breaks=brks) # y los pasamos como argumento breaks=brks
brks <- seq(from=0,to=200,by = 0.00005) # Mientras ambos histogramas tengan el mismo ancho de barra, se verán con la misma definición
hist(nodos$tau_prior,xlim=c(0,.005),breaks=brks)
# Siempre pueden consultar mínimos y máximos de los nodos "problemáticos" con nuestra
# buena amiga "summary()":
summary(nodos$tau_prior)
V1 Min. : 0.00028 1st Qu.: 0.00050 Median : 0.00114 Mean : 0.21677 3rd Qu.: 0.00479 Max. :169.45781
# Quitando las barreras negras de las barras "border=F", y colocando ambos histogramas en la misma gráfica (add=T en el segundo histograma)
brks <- seq(from=0,to=200,by = 0.00005)
hist(nodos$tau_post,xlim=c(0,.005),breaks=brks,col='#165f28',border=F)
hist(nodos$tau_prior,xlim=c(0,.005),breaks=brks,col='#31a62e',border=F,add=T)
# Extrayendo información del histograma para presentarla en otros formatos sin barras:
brks <- seq(from=0,to=200,by = 0.00005)
hist(nodos$tau_post,breaks=brks,plot=F)->ht_post
hist(nodos$tau_prior,breaks=brks,plot=F)->ht_prior
# Primer plot
# Líneas sin relleno
plot(NULL, # Abre gráfica vacía (pero para que funcione tenemos que especificar sus límites:)
# En Y, el inferior es 0 y el superior es el máximo de los máximos de densidad
# (esto garantiza que ambas distribuciones cabrán en la gráfica)
ylim=c(0,max(max(ht_prior$density),max(ht_post$density))),
# En X podemos jugar con valores arbitrarios hasta encontrar
# los que presenten la mejor comparación entre ambas distribuciones:
xlim=c(0,0.005))
lines(ht_prior$mids,ht_prior$density,col='#165f28',lwd=6)
lines(ht_post$mids,ht_post$density,col='#31a62e',lwd=6)
# Segundo plot
# Agregando polígonos de relleno
plot(NULL, # Abre gráfica vacía (pero para que funcione tenemos que especificar sus límites:)
# En Y, el inferior es 0 y el superior es el máximo de los máximos de densidad
# (esto garantiza que ambas distribuciones cabrán en la gráfica)
ylim=c(0,max(max(ht_prior$density),max(ht_post$density))),
# En X podemos jugar con valores arbitrarios hasta encontrar
# los que presenten la mejor comparación entre ambas distribuciones:
xlim=c(0,0.005))
polygon(x=c(head(ht_prior$mids,1),ht_prior$mids,tail(ht_prior$mids,1)),
y=c(0,ht_prior$density,0),
col='#165f2855')
polygon(x=c(head(ht_post$mids,1),ht_post$mids,tail(ht_post$mids,1)),
y=c(0,ht_post$density,0),
col='#31a62e55')
lines(ht_prior$mids,ht_prior$density,col='#165f28',lwd=6)
lines(ht_post$mids,ht_post$density,col='#31a62e',lwd=6)
# Coloreando barras de acuerdo con su densidad/altura relativa
# (¡utilizando una función lineal!)
brks <- seq(from=0,to=200,by = 0.00005)
hist(nodos$tau_prior,breaks=brks,plot=F)->ht
b1 <- 255/max(ht$density)
b0 <- 0
rgb_h <- b1*ht$density+b0
colores <- rgb(red = rgb_h,blue = 120,green = 100,maxColorValue = 255)
hist(nodos$tau_prior,xlim=c(0,0.005),breaks=brks,col=colores)